Historic and contemporary biogeographic perspectives on range‐wide spatial genetic structure in a widespread seagrass

Abstract Historical and contemporary processes drive spatial patterns of genetic diversity. These include climate‐driven range shifts and gene flow mediated by biogeographical influences on dispersal. Assessments that integrate these drivers are uncommon, but critical for testing biogeographic hypotheses. Here, we characterize intraspecific genetic diversity and spatial structure across the entire distribution of a temperate seagrass to test marine biogeographic concepts for southern Australia. Predictive modeling was used to contrast the current Posidonia australis distribution to its historical distribution during the Last Glacial Maximum (LGM). Spatial genetic structure was estimated for 44 sampled meadows from across the geographical range of the species using nine microsatellite loci. Historical and contemporary distributions were similar, with the exception of the Bass Strait. Genetic clustering was consistent with the three currently recognized biogeographic provinces and largely consistent with the finer‐scale IMCRA bioregions. Discrepancies were found within the Flindersian province and southwest IMCRA bioregion, while two regions of admixture coincided with transitional IMCRA bioregions. Clonal diversity was highly variable but positively associated with latitude. Genetic differentiation among meadows was significantly associated with oceanographic distance. Our approach suggests how shared seascape drivers have influenced the capacity of P. australis to effectively track sea level changes associated with natural climate cycles over millennia, and in particular, the recolonization of meadows across the Continental Shelf following the LGM. Genetic structure associated with IMCRA bioregions reflects the presence of stable biogeographic barriers, such as oceanic upwellings. This study highlights the importance of biogeography to infer the role of historical drivers in shaping extant diversity and structure.

Glacially mediated sea level changes affected the extent of continental shelf habitat available, where habitats were eliminated or significantly reduced, and reformed and expanded over time (Dolby et al., 2016). The lowest recorded sea level at ~120 m below present occurred during the Last Glacial Maximum (LGM) between 21,000 and 19,000 years BP (Lewis et al., 2013), leaving most of Australia's current continental shelf exposed, and a land bridge connecting mainland Australia to Tasmania. Shoreline reconstructions show marine inundation of the Bassian Landbridge began ~17,500 years ago from the west (figure 4 in Lambeck & Chappell, 2001). Drowned (submerged) shorelines persist (Brooke et al., 2017), leaving a paleorecord as rising sea levels inundated the shelf. Independent verification using carbon dating of P. australis sheath remains from seafloor sedimentary cores in western and central-southern Australia, estimated at up to 3500 years old (Serrano et al., 2016), is consistent with recolonization of inshore coastal environments post-LGM. However, these clues provide no insight as to whether recolonization occurred from discrete refugia or more widespread edge-of-shelf meadows.
These are based on species turnover, and a combination of extrinsic (water movement and historical barriers) and intrinsic (life history and habitat type) factors that explain why these relationships largely hold across a range of taxa (e.g., Ayre et al., 2009;Li et al., 2013;York et al., 2008), despite the absence of barriers, such as the Bassian Landbridge (figure 2 in Williams et al., 2018). These broad-scale provinces are, however, likely to miss biologically meaningful barriers that have caused significant genetic structure. In response, Waters et al. (2010) recommended a quantitative approach to testing finer-scale, regional marine biogeographical frameworks.

An Integrated Marine and Coastal Regionalisation Framework for
Australian waters (IMCRA v4.0;Commonwealth of Australia, 2006;Last et al., 2010) was developed to capture spatial patterns in species distributions. The IMCRA framework defined inshore coastal provincial bioregions based on demersal fish biogeography, of which 10 regions span southern Australia (Figure 1). Six are classified as regions of "biotic endemism" and include subtropical, warm temperate, or cool temperate waters, and four are identified as "transitions" which are less well-defined mixing areas. This regional framework can be used to classify populations and communities into bioregions that make ecological sense (e.g., reflect connectivity) and are at a scale useful for regional planning. A multispecies approach by Pope et al. (2015) determined that when latitudinal information was combined with the IMCRA model for 101 species, genetic diversity in populations increased toward the equator, although a notable "hump" with high genetic diversity was observed at the temperatetropical zone along the Western Australian coastline (IMCRA 29; Pope et al., 2015). Here, we test whether these bioregions are biologically meaningful for a widespread temperate benthic habitatforming seagrass species, Posidonia australis.
The genus Posidonia belongs to an ancient group of marine angiosperms, estimated to have diverged more than 50 mya, based on molecular fossil-dated phylogenies (Waycott et al., 2018). The extant Australian members of the genus began diverging approximately 12.8 mya (Aires et al., 2011) and possess similar biological attributes including a perennial habit capable of forming clonal plants with potential for extreme longevity and size (e.g., Arnaud-Haond et al., 2012;Edgeloe et al., 2022). Ribbon weed, Posidonia australis Hook.f., is endemic to the temperate waters of southern Australia (Edgar, 2000). Its distribution spans approximately 5300 km of coastline south from Shark Bay in the temperate-tropical zone in Western Australia to Wallis Lake in central New South Wales. Distribution records show the species is naturally disjunct across its range, and increasing anthropogenic impacts have caused significant declines (Evans et al., 2018;Short et al., 2011). This species typically grows in large continuous meadows from the low tide mark to 15 m depth, favoring more sheltered bays and estuaries (Cambridge & Kuo, 1979;Carruthers et al., 2007). Fruit production is highly variable across the species' range, but often prolific in the southwest of Australia.
The fruit, containing a single direct-developing seed with no dormancy, begins releasing from lower latitudes in the austral spring (November-January; Kuo & McComb, 1989). Positively buoyant fruit are dispersed on the sea surface via currents and windage , frequently covering distances of 10 s of km over a few days (Ruiz-Montoya et al., 2015;Sinclair et al., 2018).
We modeled the historical and contemporary distributions of P. australis and assessed range-wide genetic diversity using microsatellite DNA data to address the following questions: (1) How is genetic diversity spatially structured among meadows sampled from across the 5300 km range of this species; (2) Does this spatial genetic structure correspond to the broad-and fine-scale biogeographical models for the temperate Australian continental shelf; and (3) Were historical meadows widespread or restricted in isolated refugia during the LGM? 2 | MATERIAL S AND ME THODS

| Field sampling and genotyping
This study includes collections made from across the species' entire geographic range between 2004 and 2016. These meadows covered the three biogeographic provinces and 8 of 10 IMCRA bioregions (Table 1; Figure 1). Thirty individual shoots from each of the 44 P. australis meadows were sampled following methods described by Sinclair et al. (2014). DNA was extracted from meristem tissue and genotyped for nine polymorphic microsatellite loci (PaA1, PaA105, PaA120, PaB6, PaB8, PaB112, PaD12, PaD113, and Pa118/9), using methods previously described (Sinclair et al., 2009. We combined data from published regional studies (Evans et al., 2014;Sinclair et al., 2014Sinclair et al., , 2016 with sampling from an additional 14 meadows for a range-wide assessment. A subset of samples were used as positive controls in each run to ensure consistent scoring of alleles across locations. There was no widespread evidence for linkage disequilibrium, presence of null alleles, stuttering, or large allele dropout from our previous studies, as assessed using MICRO-CHECKER (van Oosterhout et al., 2004).
Patterns of allelic diversity were assessed visually (presence of consecutive-sized alleles) to determine whether individual loci followed a stepwise model of mutation (addition or subtraction of single repeat units; Kimura & Ohta, 1978). Linear regressions were used to test the relationship between latitude (north to south along both west and east Australian coasts) and longitude (west to east across southern Australia) with measures of genetic diversity: clonal diversity (R), allelic diversity (Na), and expected heterozygosity (H e ).
A principal coordinate analysis (PCoA) was implemented in Genalex to visualize the genetic relationships among sampled meadows based on the complete dataset and repeated on a reduced dataset containing only unique MLGs. A Bayesian clustering approach was implemented in STRUCTURE v2.3.4 (Pritchard et al., 2000) to infer significant genetic clusters based on the complete dataset only.

A reduced dataset (MLGs only) resulted in an uneven number of
MLGs within sampled meadows due to extensive clonality in some meadows, which can create biases (Puechmaille, 2016). This approach identified the number of K clusters (or populations), assigned individuals to them, and identified admixed or migrant individuals.
The analysis was performed for K = 1-44. A total of 10 independent runs were made for each value of K with 10,000 "burn-in" and 100,000 replicates, using an admixture model with allele frequencies correlated across populations and without sampling locations specified as priors. Calculation of ΔK was used to infer the most likely number of clusters for 10 replicate runs for each K (Evanno et al., 2005) implemented using STRUCTURE HARVESTER v6.94 (Earl & von Holdt, 2012). A STRUCTURE bar plot was drawn using Structure Plot v2.0 (Ramasamy et al., 2014). Population differentiation was also estimated within each significant K cluster to determine regional structure, F ST (Wright, 1943) and D (Jost, 2008), in Genalex.
A hierarchical analysis of molecular variation (AMOVA) was performed in Genalex to test for significance (at different spatial scales) among a priori recognized biogeographical provinces (Flindersian, Maugean, and Peronian), IMCRA inshore coastal provincial bioregions (29-38), and three "soft" seascape features as barriers to dispersal (Capes Upwelling, Great Australian Bight, and Bonney Upwelling; Figure 1), relative to genetic clustering. Components of variance were computed at these multiple levels, with significance based on 9999 permutations.
Overall and pairwise population differentiation was estimated using F ST (Wright, 1943) and D (Jost, 2008) in Genalex. Isolationby-distance (IBD) relationships were assessed through the association between genetic differentiation (F ST and D) and oceanographic distance (km) across the whole species range and by coastlines: latitudinal eastern (12 meadows; 32-37°S), western (12 meadows; 25-35°S), and longitudinal southern (20 meadows; 117-148°E) using Mantel tests (9999 iterations). Oceanographic distance was calculated using the "gdistance" package in R (3.4.4), which estimated the least-cost (shortest) distance via water between GPS locations of all sampled meadows.

| Mapping the historical and current distribution
Predictive modeling techniques were used to map the current distribution of P. australis, while the historical distribution of Posidonia habitat was created using depth as the primary distribution driver. Presence records for P. australis were downloaded from the Global Biodiversity Information Facility (GBIF: https://www.gbif.org/en/), of which 1153 reliable records were mapped using ArcGIS, after duplicates and points with errors in geo-positioning were removed ( Figure S1). A total of 23 environmental parameters that were potential predictors of P. australis occurrence were compiled from Bio-ORACLE (Tyberghein et al., 2012; http://www.bio-oracle.org/) and GeoScience Australia (https://www. ga.gov.au/data-pubs) online data portals, at a spatial resolution of 5 arcmin. Highly correlated variables were removed to avoid excessive autocorrelation among predictors and reduce model overfitting based on variance inflation factor (VIF) values, which was calculated for each variable using the "vifstep" function in the "usdm" R package (Naimi et al., 2014). A threshold of VIF > 5 was used in this study as a conservative measure (Naimi et al., 2014). A final list of nine variables was used to build models (Table S1). Models were generated using five different methods and a model ensemble in the "sdm" R package (Naimi & Araújo, 2016 individual models for each method. Models were evaluated using AUC (area under the curve of a receiver-operator characteristic plot).
Model ensembles were then created in the "sdm" package (Naimi & Araújo, 2016). Models with a test AUC > 0.7 were interpreted as a good performance (e.g., Coetzee et al., 2009;Swets, 1988), and were used to create the model ensemble which was then used to project P. australis probability of occurrence across temperate Australia. A conservative presence threshold of 0.7 was used to map the final distribution. The historical distribution for P. australis was predicted based on current knowledge of light requirements, bathymetry, and paleocoastline, with a potential distribution to a depth of 20 m.

| Genetic diversity within meadows
Tri-allelic genotypes were observed in some meadows (12.9% of samples; Table 1). These were reduced to diploid genotypes for all analyses as follows. Alleles that were not detected in a homozygous form, or were rare (f < 0.05), were removed. Identical, commonly occurring triallelic genotypes were reduced to the same diploid genotype, so as not to alter the total number of MLGs (as per Sinclair et al., 2020).

| Spatial genetic structure and biogeography
Overall spatial genetic structure assessed using PCoA was similar using the complete dataset and MLGs only, despite the variable number of MLGs within meadows. The PCoA grouped meadows into five genetic clusters, consistent with spatial clustering of collection records in the GBIF database (Figures 3, S1) A plot of the absolute second-order rate of change of the likelihood distribution from the STRUCTURE analysis of the complete data showed a large peak at K = 2, which we attribute to the K = 2 conundrum (Janes et al., 2017). The "optimal" number of K clusters = 5 At a higher level, genetic clustering was largely consistent with the three a priori marine biogeographic provinces (AMOVA, p < .001; Table 2). The Maugean (Bass Strait) and Peronian (east coast) genotypes were more closely related to each other than the Flindersian ( Figure 3). Genetic structuring was also largely consistent with the finer-scale IMCRA bioregions, accounting for 33% of variation (  Table 2).
Overall genetic differentiation across the species range was high

| Historical and current distribution
The historical and current distribution models predict mostly con-

| DISCUSS ION
Hierarchical patterns of range-wide spatial genetic structure in P.

| Long-term stability and broad biogeographic patterns
The combination of predictive mapping and/or population genetic assessments from multiple studies (e.g., Teske et al., 2017Teske et al., , 2018Nielsen et al., 2021;Assis et al., 2022; this study) leads us to propose several broad biogeographical hypotheses for southern hemisphere temperate marine species: (1) strong latitudinal gradient in intraspecific diversity (highest diversity in southern populations); (2) past distributions extended into lower latitudes (more northerly distribution); (3) widespread distributions of coastal shelf species persisted during the LGM (not distinct, isolated refugia, as proposed by Waters and Roy (2003)); and (4) important role of "soft" oceanic features such as currents and upwellings, as shared historical and contemporary barriers to genetic connectivity. These are in striking contrast to equivalent latitudes in the northern hemisphere where glaciation caused marine taxa to retreat to periglacial and southern refugia (Maggs et al., 2008). Zostera marina, also highlights the historical influence on modern marine ecosystems (Duffy et al., 2022).
The largely similar widespread historical and current Posidonia distributions suggest shared habitat connectivity over geological timescales (Miocene to present), with common seascape features influencing distributions and genetic structure. Two notable historical differences included range extensions into lower latitudes and the Bassian Landbridge (Figures 1 and 5). This long-term persistence in habitat connectivity is likely driven by a stable Australian continental shelf and offshore boundary currents creating an environment that favored the persistence of widespread shallow-water benthic communities. Extended periods of lower sea levels reduced the continental shelf area available for benthic marine ecosystems (Heap & Harris, 2008;James & Bone, 2011;Williams et al., 2018) and left a record of submerged paleoshorelines (Brooke et al., 2017;James & Bone, 2011). Given the high light requirement of Posidonia, the shallow, more gently sloping continental shelf margins would have likely retained large, TA B L E 2 Test of biogeographic hypotheses and three marine features: hierarchical analysis of molecular variance among sampled Posidonia australis meadows. well-connected meadows, as seen in current distributions and supported by the historical distribution model. Contemporary shelf habitats were then recolonized by widespread seagrass meadows tracking rising sea levels. For example, genetic data for P. australis ( Figure 2; Sinclair et al., 2016) suggest initial recolonization of the Bass Strait from the upper slope of the Otway Shelf to the west, which was also found for a sea star, Coscinasterias muricata (Waters & Roy, 2003). The secondary contact or admixture zones between southern and east coast meadows are evidence of this biogeographical barrier, which has been present on multiple occasions through millennia (>67 m below current sea level ;Lambeck & Chappell, 2001).

Source of variation
In striking contrast to the general finding of an inverse relationship between genetic diversity and latitude for Australian marine species (Pope et al., 2015), we found that genetic diversity increased Spatial patterns of genetic structure across the world's longest boundary circulation, which traverses Australia's southern coastline (Middleton & Bye, 2007), were weakly associated with two habitat-forming species, P. australis (this study) and the kelp Ecklonia radiata (Coleman et al., 2011). Ecklonia radiata appeared to be correlated with the strength of boundary currents (Coleman et al., 2011), while the pattern in P. australis was more complex.
Sparse sampling across southern Australia has often limited wider conclusions (Coleman et al., 2011;Kassahn et al., 2003;Stiller et al., 2021;Wilson et al., 2017), however, the higher genetic diversity in P. australis meadows may be associated with a change in ploidy, as recently detected in Shark Bay meadows (Edgeloe et al., 2022).

| Provincial bioregions and spatial genetic structure
Range-wide spatial genetic structure in P. australis was largely as- Australian coast, creating a physical and temporal barrier to gene flow (see Hanson et al., 2005). This seasonal upwelling occurs during the Austral Spring/Summer, coinciding with fruit release for reproductive meadows. Significant genetic structure was also reported in this region for E. radiata (Vranken et al., 2021).
A second region of discrepancy occurred in southeastern

F I G U R E 5
The modeled current (dark green) and historical distribution of Posidonia at the Last Glacial Maximum (c. 120 m below present sea level, light green) approximates the location of the Continental shelf margin. Inset A. Larger paleomeadows were present along the margins of the Otway Shelf.
Phylogenetic structuring in leafy seadragons, Phycodurus eques (Stiller et al., 2021), also supports this relationship by suggesting that western regions were recolonized from the east during post-LGM inundation of the Continental shelf.
Posidonia australis meadows are highly fragmented along the exposed coastlines within the Southwest and Bass Strait transitions (IMCRA 30,34). The western Bass Strait transition (IMCRA 34) contains a seasonal upwelling, Bonney Upwelling (Kämpf et al., 2004) which also appears to restrict gene flow between central-southern Australia (IMCRA 33) and Bass Strait (IMCRA 35). Admixture among P. australis meadows across the southeastern transition (IMCRA 37) was reported in species with low and high dispersal capabilities (e.g., Miller et al., 2013;Stoessel et al., 2020;Waters & Roy, 2003), so there is no consensus on barriers to gene flow in this region. For example, multiple features occur including 90-mile beach, a significant change in geomorphic structure (figure 9 in Heap & Harris, 2008), and the winter 14°C sea surface isotherm (O'Hara & Poore, 2000).
The spatial genetic structure reflects the distance between available habitats, as per Ayre et al. (2009), as well as the complex coastlines and life-history traits which can create an appearance of localized "chaotic genetic patchiness" due to genetic drift and temporal variation in collective dispersal events (Broquet et al., 2013;Johnson & Black, 1982;Sinclair et al., 2014). Some seagrass species are particularly well adapted for inshore sea surface dispersal of floating fruits and/or reproductive plant parts (e.g., Hernawan et al., 2017;Sinclair et al., 2016). Contemporary dispersal by P. australis fruits is influenced by daily shifts in local prevailing winds, with distances of 100-150 km representing the upper limit for dispersing Sinclair et al., 2018). The range-wide association between genetic and geographic distance in P. australis is congruent with a widespread species displaying a regional dispersal capacity, the scale of which is broadly consistent with spatial genetic patterns found in other widespread seagrasses (reviewed in Kendrick et al., 2012).

| Implications for the future of climate change
Overall, our combination of predictive distribution mapping and analysis of spatial genetic structure has allowed us to infer how historical and contemporary seascape features have influenced the capacity of P. australis to effectively track sea level changes associated with natural climate cycles over millennia via clonal and sexual reproduction. The similarities between historical and current distributions were consistent with shared biogeographical features and long-term resilience in P. australis. However, recent increases in ocean temperatures and acidification associated with global warming create a more stressful environment, with many temperate species already on the move, retreating to deeper waters, higher latitudes, or becoming locally extinct (Pecl et al., 2017).
There is potential for a southward range expansion of P. australis down Tasmania's east coast into suitable habitat, as identified in our modeling. Projections for P. oceanica in the Mediterranean suggest it will no longer be able to inhabit its' current range, with a 75% loss of suitable habitat by 2050 and functional extinction by 2100 (Chefaoui et al., 2017). Similarly, the physiological predictions for survival of P. australis to increasing temperatures along the west coast of Australia suggest a range contraction of between 200 and 400 km by 2100 (Hyndes et al., 2016). Whole-genome duplication through polyploidy in P. australis has apparently increased thermal and salinity tolerance in extreme environments and enabled colonization of Shark Bay waters post-LGM (Edgeloe et al., 2022). This is an example of the evolutionary significance of polyploidy and shows how species and populations have and will continue to evolve and adapt to stressful environments ( Van de Peer et al., 2021). Further research is required on the role of hybridization and Holocene range expansions into the warm, hypersaline environments in the upper gulf region of central-southern Australia. Our current study shows a historic resilience to climate change in P. australis and suggests stability in the current IMCRA bioregions for interpreting biogeographic patterns. However, more detailed sampling across a range of taxa will be required to resolve many of the discrepancies, in which a nested mesoscale framework may better capture diversity and spatial genetic structure and be a preferred framework for managing biodiverse ecosystem structure.